# Analyze SSI Sample Data
# Created 01.03.2019

library(sandwich)
library(lmtest)
library(ri)
library(ggplot2)
library(RColorBrewer)
library(gridExtra)
library(xtable)
library(texreg)
library(dplyr)

rm(list = ls())

source(".../04 SSI Clean Data.R")
source(".../Go-to Stats Functions.R")

figdir <- "PATH TO FOLDER TO OUTPUT FIGURES/TABLES"

######################################
# Functions to Make Figures and Tables
######################################

## Treatment Group Means Plot Function

histmeans <- function(data, indicators, treatment, response, 
                      xlabs, ylab, title, addc, ymin, ymax) {
  
  tmeans<- as.data.frame(matrix(indicators, nrow = length(indicators)))
  names(tmeans) <- "treat"
  tmeans$treat <- factor(tmeans$treat, levels = indicators)
  
  tmeans$resp <- NA
  tmeans$cilb <- NA
  tmeans$ciub <- NA
  tmeans$pval <- NA
  tmeans$n <- NA
  
  attach(data)
  for (i in 1:length(indicators)) {
    tmeans[i, 2:6] <- tci1(data[treatment == tmeans[i, 1], response], 0.05, FALSE)
  }
  
  detach(data)
  
  tmeans$treatlab <- factor(xlabs)
  tmeans$treatlab <- factor(tmeans$treatlab, levels = xlabs)
  
  p <-  ggplot(tmeans, aes(x = treatlab, y = resp, fill = treatlab)) + 
    geom_bar(position = position_dodge(), stat ="identity") +
    geom_errorbar(aes(ymin = cilb, ymax = ciub),
                  width=.2,                    
                  position=position_dodge(.9), colour = "black") +
    
    guides(fill = FALSE) +
    xlab("") + ylab(ylab) + ggtitle(title) + coord_cartesian(ylim = c(ymin, ymax)) +
    
    scale_fill_grey(start = 0.30, end = 0.90) + 
    
    theme(panel.background = element_blank(),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          axis.text.x = element_text(size = 20, colour = "black"), axis.ticks.x = element_blank(),
          axis.text.y = element_text(size = 20, colour = "black"),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5), 
          axis.title.y = element_text(size = 20))
  
  for (j in 1:length(indicators)) {
    p <- p + annotate("text", x = j, y = tmeans[j,4] + addc, size = 7,
                      label = paste(round(tmeans[j,2], digits = 2), sep = ""))
  }
  
  return(p)
  #return(tmeans)
}

## Nutrition/Job Training Treatment Group Differences Table Function

difftable <- function(data, treatment, g) {
  
  attach(data)
  
  sampcovs <- c("% Take-up", "Cents/$ Divert Total", 
                "Prog Goal", "Admin", "Govt Misuse",
                "Citizen Misuse", 
                "less $40k v.$40-75k", "$40-75k v. $75k more", "n")
  
  resp <- c("_tu_4", "_divert", "_w_1", "_w_2", "_w_3", "_w_4", "_dist_1", "_dist_2")
  
  sampdif <- as.data.frame(matrix(NA, ncol = 3, nrow = length(sampcovs)))
  sampdif$Outcomes <- sampcovs
  
  for(i in 1:length(resp)) {
    sampdif[i, 1:3] <- tci2(data[treatment == "direct", paste(g, resp[i], sep = "")], 
                            data[treatment == "indirect", paste(g, resp[i], sep = "")], 
                            0.05, 3)[c(1,2,4)]
  }
  
  sampdif[9, 1:2] <- table(treatment)
  
  detach(data)
  
  sampdif <- sampdif[, c(4,1:3)]
  names(sampdif)[2:4] <- c("Direct", "Indirect", "p-value")
  
  # Adjust p-values
  
  sampdif[1:2, 4] <- p.adjust(sampdif[1:2, 4], "BH")
  sampdif[3:6, 4] <- p.adjust(sampdif[3:6, 4], "BH")
  sampdif[7:8, 4] <- p.adjust(sampdif[7:8, 4], "BH")
  
  return(print(xtable(sampdif), include.rownames = FALSE))
  #return(sampdif)
}

## HMID/UI Treatment Group Differences Table Function

difftable2 <- function(data, treatment, g) {
  
  attach(data)
  
  sampcovs <- c("Spillover", "Favor Status Quo", "Change Benefit", "Change Premiums",
                "% Take-up", "Cents/$ Divert Total", "Uncertainty",
                "Prog Goal", "Admin", "Govt Misuse", "Citizen Misuse", 
                "Ineq", "n")
  
  resp <- c("_r2", "_s1", "_s2", "_s3", "_tu_1", "_divert", "_uc",
            "_w_1", "_w_2", "_w_3", "_w_4", "_ineq")
  
  sampdif <- as.data.frame(matrix(NA, ncol = 3, nrow = length(sampcovs)))
  sampdif$Outcomes <- sampcovs
  
  for(i in 1:length(resp)) {
    
    if (i == 2) {
      pests <- prop.test(c(sum(data[treatment == "reg", paste(g, resp[i], sep = "")] == 3), 
                           sum(data[treatment == "pro", paste(g, resp[i], sep = "")] == 3)),
                         c(sum(data[treatment == "reg", paste(g, resp[i], sep = "")]), 
                           sum(data[treatment == "pro", paste(g, resp[i], sep = "")])),
                         correct = T)
      
      sampdif[i, 1:3] <- c(pests$estimate[1], pests$estimate[2], pests$p.value)
    }
    
    else {
      sampdif[i, 1:3] <- tci2(data[treatment == "reg", paste(g, resp[i], sep = "")], 
                              data[treatment == "pro", paste(g, resp[i], sep = "")], 
                              0.05, 3)[c(1,2,4)]
    }
  }
  
  detach(data)
  
  sampdif <- sampdif[, c(4,1:3)]
  names(sampdif)[2:4] <- c("Regressive", "Progressive", "p-value")
  
  # Adjust p-values
  
  if (g == "hm") {
    sampdif <- sampdif[-4, ]
    sampdif[2:3, 4] <- p.adjust(sampdif[2:3, 4], "BH")
    sampdif[4:6, 4] <- p.adjust(sampdif[4:6, 4], "BH")
    sampdif[7:10, 4] <- p.adjust(sampdif[7:10, 4], "BH")
    sampdif[12, 2:3] <- table(treatment)
    sampdif <- sampdif[c(2,3,11,12), ]
  }
  
  if (g == "ui") {
    sampdif[2:4, 4] <- p.adjust(sampdif[2:4, 4], "BH")
    sampdif[5:7, 4] <- p.adjust(sampdif[5:7, 4], "BH")
    sampdif[8:11, 4] <- p.adjust(sampdif[8:11, 4], "BH")
    sampdif[13, 2:3] <- table(treatment)
    sampdif <- sampdif[c(2,3,4,12,13), ]
  }
  
  return(print(xtable(sampdif), include.rownames = FALSE))
  #return(sampdif)
}

# Balance on covariates table function

pdif <- function(resp, treat, t1val, t2val) {
  
  dif <- prop.test(c(sum(resp[treat == t1val]), sum(resp[treat == t2val])), 
                   c(sum(treat == t1val), sum(treat == t2val)), correct = TRUE)
  
  result <- c(dif$estimate, dif$estimate[1] - dif$estimate[2], dif$p.value, length(treat))
  
  return(result)
}

covbal <- function(data, treat, t1val, t2val, t1lab, t2lab) {
  
  attach(data)
  
  sampcovs <- c("Age (years)", "Black (1 = yes)", "Latino (1 = yes)", "Democrat (1 = yes)", "Republican (1 = yes)", 
                "College (1 = yes)", "Income (see notes)", "n / F test")
  
  sampdif <- as.data.frame(matrix(NA, ncol = 3, nrow = length(sampcovs)))
  sampdif$Covariate <- sampcovs
  
  sampdif[sampdif$Covariate == "Age (years)", 1:3] <- tci2(age[!is.na(age) & treat == t1val], 
                                                           age[!is.na(age) & treat == t2val], 0.05, 3)[c(1,2,4)]
  
  sampdif[sampdif$Covariate == "Black (1 = yes)", 1:3] <- pdif(black[!is.na(black)], 
                                                               treat[!is.na(black)], t1val, t2val)[c(1,2,4)]
  
  sampdif[sampdif$Covariate == "Latino (1 = yes)", 1:3] <- pdif(latino[!is.na(latino)], 
                                                                treat[!is.na(latino)], t1val, t2val)[c(1,2,4)]
  
  # Dem/Rep Including Leaners
  
  sampdif[sampdif$Covariate == "Democrat (1 = yes)", 1:3] <- pdif(dem_plus_lean, treat, t1val, t2val)[c(1,2,4)]
  
  sampdif[sampdif$Covariate == "Republican (1 = yes)", 1:3] <- pdif(rep_plus_lean, treat, t1val, t2val)[c(1,2,4)]
  
  sampdif[sampdif$Covariate == "College (1 = yes)", 1:3] <- pdif(college, treat, t1val, t2val)[c(1,2,4)]
  
  sampdif[sampdif$Covariate == "Income (see notes)", 1:3] <- tci2(hhinc[treat == t1val], 
                                                                  hhinc[treat == t2val], 0.05, 3)[c(1,2,4)]
  
  sampdif[sampdif$Covariate == "n / F test", 1:2] <- table(treat)
  
  data$mod_treat <- 0
  data$mod_treat[treat == t2val] <- 1
  
  joint <- summary(lm(mod_treat ~ age + black + latino + college + dem_plus_lean + rep_plus_lean + as.factor(hhinc), data = data))
  fs <- joint$fstatistic
  fp <- pf(fs[1], fs[2], fs[3], lower.tail = FALSE)
  
  sampdif[sampdif$Covariate == "n / F test", 3] <- fp
  
  detach(data)
  
  sampdif <- sampdif[, c(4,1:3)]
  names(sampdif)[2:4] <- c(t1lab, t2lab, "p-value")
  
  return(print(xtable(sampdif), include.rownames = FALSE))
}

########################
# MAIN RESULTS FOR PAPER
########################

# Nutrition/Job Training Results

n_sup_paper <- histmeans(y, c("direct", "indirect"), n_treat_v, "n_r", 
                   c("Direct", "Indirect"), 
                   "[-2,2] = Strongly oppose...Strongly support\n", 
                   "", 0.25, 0, 1.25)
n_sup_paper

mean(y$n_r[y$n_treat_v == "indirect"]) - mean(y$n_r[y$n_treat_v == "direct"])

t.test(y$n_r[y$n_treat_v == "indirect"], y$n_r[y$n_treat_v == "direct"])

mean(y$n_r[y$n_treat_v == "indirect"]) - mean(y$n_r[y$n_treat_v == "direct"]) / 
  sd(y$n_r[y$n_treat_v == "direct"])

j_sup_paper <- histmeans(y, c("direct", "indirect"), j_treat_v, "j_r", 
                   c("Direct", "Indirect"), 
                   "[-2,2] = Strongly oppose...Strongly support\n", 
                   "", 0.25, 0, 1.25)
j_sup_paper

mean(y$j_r[y$j_treat_v == "indirect"]) - mean(y$j_r[y$j_treat_v == "direct"])

t.test(y$j_r[y$j_treat_v == "indirect"], y$j_r[y$j_treat_v == "direct"])

mean(y$j_r[y$j_treat_v == "indirect"]) - mean(y$j_r[y$j_treat_v == "direct"]) /
  sd(y$j_r[y$j_treat_v == "direct"])
  
# N/J Tables

difftable(y, n_treat_v, "n")
difftable(y, j_treat_v, "j")

# Hand Check Means: Nutrition

t.test(y$n_tu_4[y$n_treat_v == "direct"], y$n_tu_4[y$n_treat_v == "indirect"])
t.test(y$n_divert[y$n_treat_v == "direct"], y$n_divert[y$n_treat_v == "indirect"])
t.test(y$n_w_1[y$n_treat_v == "direct"], y$n_w_1[y$n_treat_v == "indirect"])
t.test(y$n_w_2[y$n_treat_v == "direct"], y$n_w_2[y$n_treat_v == "indirect"])
t.test(y$n_w_3[y$n_treat_v == "direct"], y$n_w_3[y$n_treat_v == "indirect"])
t.test(y$n_w_4[y$n_treat_v == "direct"], y$n_w_4[y$n_treat_v == "indirect"])
t.test(y$n_dist_1[y$n_treat_v == "direct"], y$n_dist_1[y$n_treat_v == "indirect"])
t.test(y$n_dist_2[y$n_treat_v == "direct"], y$n_dist_2[y$n_treat_v == "indirect"])

# Hand Check Means: Job Training

t.test(y$j_tu_4[y$j_treat_v == "direct"], y$j_tu_4[y$j_treat_v == "indirect"])
t.test(y$j_divert[y$j_treat_v == "direct"], y$j_divert[y$j_treat_v == "indirect"])
t.test(y$j_w_1[y$j_treat_v == "direct"], y$j_w_1[y$j_treat_v == "indirect"])
t.test(y$j_w_2[y$j_treat_v == "direct"], y$j_w_2[y$j_treat_v == "indirect"])
t.test(y$j_w_3[y$j_treat_v == "direct"], y$j_w_3[y$j_treat_v == "indirect"])
t.test(y$j_w_4[y$j_treat_v == "direct"], y$j_w_4[y$j_treat_v == "indirect"])
t.test(y$j_dist_1[y$j_treat_v == "direct"], y$j_dist_1[y$j_treat_v == "indirect"])
t.test(y$j_dist_2[y$j_treat_v == "direct"], y$j_dist_2[y$j_treat_v == "indirect"])

# Check robustness by subsetting to respondents who saw each policy area first

table(y$DO.BR.FL_33)
table(y$DO.BR.FL_20)
table(y$first)

table(y$n_treat_v)
table(y$j_treat_v)

table(y$first, y$n_treat_v)
table(y$first, y$j_treat_v)

table(y$DO.BR.FL_33, y$n_treat_v)
table(y$DO.BR.FL_33, y$j_treat_v)
table(y$DO.BR.FL_20, y$n_treat_v)
table(y$DO.BR.FL_20, y$j_treat_v)

n_sup_paper_rob <- histmeans(y[y$DO.BR.FL_33 == "Nutrition|Job Training" | y$DO.BR.FL_20 == "Nutrition|Job Training", ],
                             c("direct", "indirect"), n_treat_v, "n_r", 
                         c("Direct", "Indirect"), 
                         "[-2,2] = Strongly oppose...Strongly support\n", 
                         "", 0.25, 0, 1.25)


j_sup_paper_rob <- histmeans(y[y$DO.BR.FL_33 == "Job Training|Nutrition" | y$DO.BR.FL_20 == "Job Training|Nutrition", ],
                             c("direct", "indirect"), j_treat_v, "j_r", 
                         c("Direct", "Indirect"), 
                         "[-2,2] = Strongly oppose...Strongly support\n", 
                         "", 0.25, 0, 1.25)

n_sup_paper_rob
j_sup_paper_rob

difftable(y[y$DO.BR.FL_33 == "Nutrition|Job Training" | y$DO.BR.FL_20 == "Nutrition|Job Training", ], n_treat_v, "n")
difftable(y[y$DO.BR.FL_33 == "Job Training|Nutrition" | y$DO.BR.FL_20 == "Job Training|Nutrition", ], j_treat_v, "j")

# HMID/UI Results

hm_r1_dist <- histmeans(y[y$policy == "hm" & !is.na(y$hm_treat_d), ], c("reg", "pro"), hm_treat_d, 
                        "hm_r1", c("Regressive", "Progressive"), 
                       "Dec a lot = -1, Same = 0, Inc a lot = 1\n",
                       "", 0.1, 0, 0.4)

hm_r1_dist

t.test(y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "reg", "hm_r1"], 
       y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "pro", "hm_r1"])

###
robreg(lm(hm_r1 ~ hm_treat_v, data = y[y$policy == "hm",]))
robreg(lm(hm_r1 ~ hm_treat_d + hm_treat_v + hm_treat_d*hm_treat_v, data = y[y$policy == "hm" & !is.na(y$hm_treat_d),]))
###

ui_r1_dist <- histmeans(y[y$policy == "ui" & !is.na(y$ui_treat_d), ], c("reg", "pro"), ui_treat_d, 
                        "ui_r1", c("Regressive", "Progressive"), 
                        "Dec a lot = -1, Same = 0, Inc a lot = 1\n",
                        "", 0.1, 0, 0.4)

ui_r1_dist

t.test(y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "reg", "ui_r1"], 
       y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "pro", "ui_r1"])

hm_pc1_dist <- histmeans(y[y$policy == "hm" & !is.na(y$hm_treat_d), ], c("reg", "pro"), hm_treat_d,
                         "hm_s2", c("Regressive", "Progressive"), 
                        "[-1,1] = Low Less/High More...Low More/High Less\n", 
                        "", 0.25, 0, 1)
hm_pc1_dist

ui_pc1_dist <- histmeans(y[y$policy == "ui" & !is.na(y$ui_treat_d), ], c("reg", "pro"), ui_treat_d,
                         "ui_s2", c("Regressive", "Progressive"), 
                        "[-1,1] = Low More/High Less...Low Less/High More\n", 
                        "", 0.25, 0, 1)
ui_pc1_dist

ui_pc2_dist <- histmeans(y[y$policy == "ui" & !is.na(y$ui_treat_d), ], c("reg", "pro"), ui_treat_d,
                         "ui_s3", c("Regressive", "Progressive"),
                        "[-1,1] = Low Less/High More...Low More/High Less\n", 
                        "", 0.25, 0, 1)
ui_pc2_dist

# HMID/UI Regressive/Progressive Tables

difftable2(y[y$policy == "hm" & !is.na(y$hm_treat_d), ], hm_treat_d, "hm")
difftable2(y[y$policy == "ui" & !is.na(y$ui_treat_d), ], ui_treat_d, "ui")

# Hand Check Means/Proportions: HM & UI
# "_s1", "_s2", "_s3", "_ineq"

prop.test(c(sum(y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "reg", "hm_s1"] == 3), 
            sum(y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "pro", "hm_s1"] == 3)),
          c(sum(y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "reg", "hm_s1"]), 
            sum(y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "pro", "hm_s1"])),
          correct = T)

prop.test(c(sum(y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "reg", "ui_s1"] == 3), 
            sum(y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "pro", "ui_s1"] == 3)),
          c(sum(y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "reg", "ui_s1"]), 
            sum(y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "pro", "ui_s1"])),
          correct = T)

t.test(y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "reg", "hm_s2"], 
       y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "pro", "hm_s2"])

t.test(y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "reg", "hm_ineq"], 
       y[y$policy == "hm" & !is.na(y$hm_treat_d) & y$hm_treat_d == "pro", "hm_ineq"])

t.test(y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "reg", "ui_s2"], 
       y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "pro", "ui_s2"])

t.test(y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "reg", "ui_s3"], 
       y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "pro", "ui_s3"])

t.test(y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "reg", "ui_ineq"], 
       y[y$policy == "ui" & !is.na(y$ui_treat_d) & y$ui_treat_d == "pro", "ui_ineq"])

# MTurk Pilot Data

source(".../00 MTurk Clean (run via 05 script).R")

# Data is read in and stored as "x"

names(x)
table(x$treatd, useNA = "always")
table(x$treatv, useNA = "always")

sampcovs <- c("Support", "Spillover", "n")
sampdif <- as.data.frame(matrix(NA, ncol = 4, nrow = length(sampcovs)))
sampdif$Outcomes <- sampcovs

r1_mt_reg <- robreg(lm(r1 ~ as.factor(treatv) + as.factor(treatd),
                    subset = policy == "hm", data = x))

r1_mt_reg2 <- robreg(lm(r1 ~ as.factor(treatv),
                       subset = policy == "hm" & treatd == "reg", data = x))

r4_mt_reg <- robreg(lm(r4 ~ as.factor(treatv) + as.factor(treatd),
                    subset = policy == "hm", data = x))

sampdif[1, 3:4] <- c(r1_mt_reg$coefficients[2], r1_mt_reg$rob.pvals[2])
sampdif[2, 3:4] <- c(r4_mt_reg$coefficients[2], r4_mt_reg$rob.pvals[2])

r1_ssi_reg <- robreg(lm(hm_r1 ~ as.factor(hm_treat_v),
                       subset = policy == "hm" & !is.na(hm_treat_v), data = y))

r2_ssi_reg <- robreg(lm(hm_r2 ~ as.factor(hm_treat_v),
                        subset = policy == "hm" & !is.na(hm_treat_v), data = y))

sampdif[1, 1:2] <- c(r1_ssi_reg$coefficients[2], r1_ssi_reg$rob.pvals[2])
sampdif[2, 1:2] <- c(r2_ssi_reg$coefficients[2], r2_ssi_reg$rob.pvals[2])

sampdif <- sampdif[, c(5,1:4)]
names(sampdif)[2:5] <- c("ATE-SSI", "p-value", "ATE-MTurk", "p-value")

sampdif[3, c(2, 4)] <- c(nrow(y[y$policy == "hm" & !is.na(y$hm_treat_v), ]),
                         sum(x$policy == "hm"))

print(xtable(sampdif), include.rownames = FALSE)

### ATE-MTurk as % SD of Direct

r1_mt_reg$coefficients[2]/sd(x$r1[x$policy == "hm" & x$treatv == "direct"])

#################################
# Supplemental Information Tables
#################################

# data, treat, t1val, t2val, t1lab, t2lab

covbal(y, n_treat_v, "direct", "indirect", "Direct", "Indirect")
covbal(y, j_treat_v, "direct", "indirect", "Direct", "Indirect")

covbal(y[!is.na(y$hm_treat_d) & y$policy == "hm", ], hm_treat_d, "reg", "pro", "Regressive", "Progressive")
covbal(y[!is.na(y$ui_treat_d) & y$policy == "ui", ], ui_treat_d, "reg", "pro", "Regressive", "Progressive")

covbal(y[!is.na(y$hm_treat_v) & y$policy == "hm", ], hm_treat_v, "direct", "indirect", "Direct", "Indirect")

# Run Table for MTurk Sample

x$dem_plus_lean <- x$dem
x$rep_plus_lean <- x$rep

table(x$ps3, useNA = "always")
x$hhinc <- NA
x$hhinc[!is.na(x$ps3) & x$ps3 == 1] <- 10
x$hhinc[!is.na(x$ps3) & x$ps3 == 2] <- 30
x$hhinc[!is.na(x$ps3) & x$ps3 == 3] <- 57.5
x$hhinc[!is.na(x$ps3) & x$ps3 == 4] <- 87.5
x$hhinc[!is.na(x$ps3) & x$ps3 == 5] <- 150
x$hhinc[!is.na(x$ps3) & x$ps3 == 6] <- 200

table(x$race, useNA = "always")
x$black <- 0
x$black[!is.na(x$race) & x$race == 2] <- 1
table(x$race, x$black, useNA = "always")

table(x$race, useNA = "always")
x$latino <- 0
x$latino[!is.na(x$race) & x$race == 3] <- 1
table(x$race, x$latino, useNA = "always")

table(x$edu, useNA = "always")
x$college <- 0
x$college[x$edu > 4] <- 1
table(x$edu, x$college)

table(x$born, useNA = "always")
x$age <- NA
x$age <- 2017 - (x$born + 1919)
table(x$born, x$age)

table(x$sex, useNA = "always")
x$female <- NA
x$female[!is.na(x$sex) & x$sex == 1] <- 0
x$female[!is.na(x$sex) & x$sex == 2] <- 1
table(x$sex, x$female, useNA = "always")

covbal(x[x$policy == "hm", ], treatv, "direct", "indirect", "Direct", "Indirect")

# Sample Characteristics
# SSI

names(y)

covs <- c("age", "female", "black", "latino", "dem_plus_lean", "rep_plus_lean", "college", "hhinc")
covs_labs <- c("Age", "Female (Yes = 1)", "Black (Yes = 1)", "Latino (Yes = 1)", 
               "Democrats (Yes = 1)", "Republicans (Yes = 1)",
               "College (Yes = 1)", "Income (see notes)")

dssi <- data.frame(Covariate = covs_labs, var = covs, 
                 Mean = rep(NA, length(covs)), SD = rep(NA, length(covs)))

for (i in 1:length(covs)) {
  dssi[i,3] <- mean(y[ ,covs[i]], na.rm = T)
  
  if (dssi[i,2] %in% c("age", "hhinc")) {
    dssi[i,4] <- sd(y[ ,covs[i]], na.rm = T)
  }
}

print(xtable(select(dssi, -var)), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ssi_demos", ".tex", sep = ""))

names(x)

dmturk <- data.frame(Covariate = covs_labs, var = covs, 
                   Mean = rep(NA, length(covs)), SD = rep(NA, length(covs)))

for (i in 1:length(covs)) {
  dmturk[i,3] <- mean(x[ ,covs[i]], na.rm = T)
  
  if (dmturk[i,2] %in% c("age", "hhinc")) {
    dmturk[i,4] <- sd(x[ ,covs[i]], na.rm = T)
  }
}

print(xtable(select(dmturk, -var)), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "mturk_demos", ".tex", sep = ""))
